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We report the first experimental test of an analytic image reconstruction algorithm for optical 
tomography with large data sets. Using a continuous- wave optical tomography system with 10 8 
source- detector pairs, we demonstrate the reconstruction of an absorption image of a phantom 
consisting of a highly-scattering medium with absorbing inhomogeneities. 



Optical tomography (OT) is a biomedical imaging 
modality that utilizes diffuse light as a probe of tissue 
structure and function 1]. Clinical applications include 
imaging of breast disease and functional neuroimaging. 
The physical problem that is considered is to reconstruct 
the optical properties of an inhomogenous medium from 
measurements taken on its surface. In a typical experi- 
ment, optical fibers are used for illumination and detec- 
tion of the transmitted light ^ S 3 The number of 
measurements (source-detector pairs) which can be ob- 
tained, in practice, varies between 10 2 — 10 4 . A recently 
proposed alternative to fiber-based experiments is to em- 
ploy a narrow incident beam for illumination. The beam 
can be scanned over the surface of the medium while a 
lens-coupled CCD detects the transmitted light. Using 
such a "noncontact" method, it is possible to avoid many 
of the technical difficulties which arise due to fiber-sample 
interactions 0, 0, 0, Q • I n addition, extremely large data 
sets of approximately 10 8 — 10 10 measurements can read- 
ily be obtained. Data sets of this size have the potential 
to vastly improve the quality of reconstructed images in 
OT. 

The reconstruction of images from large data sets is 
an extremely challenging problem due to the high com- 
putational complexity of numerical approaches to the in- 
verse problem in OT. To address this challenge, we have 
developed analytic methods to solve the inverse prob- 
lem |j§ [ll|, [l2| . These methods lead to a dramatic reduc- 
tion in computational complexity and have been applied 
in numerical simulations to data sets as large as 10 10 mea- 
surements [llj. In this Letter, we report the first experi- 
mental test of an analytic image reconstruction method. 
By employing a noncontact OT system with 10 8 source- 
detector pairs, we reconstruct the optical absorption of a 
highly-scattering medium. The results demonstrate the 
feasibility of image reconstruction for OT with large data 
sets. 

We begin by considering the propagation of diffuse 
light. The density of electromagnetic energy u(r) in an 
absorbing medium obeys the diffusion equation 



DV 2 u{y) + a(r)u(r) = S(r) 



(1) 



sion constant. The energy density also obeys the bound- 
ary condition u + £h • Vu = on the surface bounding 
the medium, where n is the unit outward normal and i 
is the extrapolation length |l0|. The relative intensity 
measured by a point detector at r 2 due to a point source 
at ri is given, within the accuracy of the first Rytov ap- 
proximation, by the integral equation 



0(n,r 2 ) = J d 3 rG(r u r)G(r,r 2 )Sa(r) , 



(2) 



where the source and detector are oriented in the in- 
ward and outward normal directions, respectively [T^ . 
Here Sa(r) = a(r) — ao denotes the spatial fluctua- 
tions in a(r) relative to a reference medium with ab- 
sorption ao, G is the Green's function for Eq. (QJ with 
a = ao, and the data function <\> is defined by </>(i*i, r 2 ) = 
-G(ri,r 2 )ln(/(ri,r2)//o(ri,r 2 )), where J(ri,r 2 ) de- 
notes the intensity in the medium and io(ri,r 2 ) is the 
intensity in the reference medium. Note that the inten- 
sity is related to the Green's function by the expression 



/(ri ' r2) = 17 1 + 7 



G(n,r 2 ) 



(3) 



where a(r) is the absorption coefficient, S(r) is the power 
density of a continuous wave source, and D is the diffu- 



where So is the source power and the transport mean 
free path £* is related to the diffusion coefficient by D = 
l/3cT. 

We have constructed a noncontact OT system to 
test the analytic method of image reconstruction. 
A schematic of the instrument is shown in Fig. 1. 
The source is a continuous-wave stabilized diode laser 
(DL7140-201, Thorlabs) operating at a wavelength of 785 
nm with an output power of 70 mW. The laser output 
is divided into two beams by a beam splitter. The re- 
flected beam is incident on a power meter which moni- 
tors the stability of the laser intensity. The transmitted 
beam passes through a lens onto a pair of galvanometer- 
controlled mirrors (SCA 750, Lasesys). The mirrors are 
used to scan the beam, which has a focal spot size of 
200 /im, in a raster fashion over the surface of the sam- 
ple. After propagating through the sample, the trans- 
mitted light passes through a band-pass interference filter 
(10LF20-780, Newport) and is imaged onto a front illumi- 
nated thermoelectric-cooled 16-bit CCD array (DV435, 
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Andor Technology) using a 23 mm//1.4 lens. A mechan- 
ical shutter is placed in front of the CCD to reduce arti- 
facts associated with frame transfer within the CCD chip. 
A pulse generator with digital delay is used to trigger and 
synchronize the CCD, the shutter and the position of the 
beam. 
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Fig. 1. 

system. 



Schematic of the noncontact optical tomography 



The sample chamber is a rectangular box of depth 5 
cm with square faces of area 50 x 50 cm 2 constructed of 
clear acrylic sheets. The beam is scanned on one face of 
the sample and the opposite face is imaged by the CCD. 
The chamber is placed equidistantly from the CCD and 
the laser source along the optical axis at a distance of 
110 cm. The chamber is filled with a scattering medium 
which consists of a suspension of 1% Intralipid in water 
in which absorbing objects may be suspended. 

A tomographic data set is acquired by raster scanning 
the beam over a 29 x 29 square lattice with a lattice 
spacing of 0.5 cm. This yields 841 source positions within 
a 14 x 14 cm 2 area centered on the optical axis. For each 
source, a 429 x 429 pixel region of interest is read out from 
the CCD. This results in 184,041 detectors arranged in a 
square lattice with an effective lattice spacing equivalent 
to 0.065 cm and all detectors located within a 28 x 28 
cm 2 area centered on the optical axis. Thus a data set 
of 1.5 x 10 8 source-detector pairs is acquired. 

The inverse problem in OT consists of reconstructing 
5a from measurements of <fi. In this Letter, we consider 
the inversion of the integral equation © in the slab mea- 
surement geometry. The approach taken is to construct 
the singular value decomposition of the integral operator 
whose kernel is defined by @ and to use this result to 
obtain the pseudoinverse solution to (0). 

The starting point for this development is to consider 
the lattice Fourier transform of the sampled data function 
which is defined by 

0(qi,q2) = ^2 exp t^ qi ' ri + q2 ' r2 )l ^( ri ' r2 ) ' ( 4 ) 

ri,r 2 

where the sum is carried out over the square lattices of 
sources and detectors with lattice spacings hi and re- 
spectively. The wave vectors qi and q2 belong to the 



first Brillouin zones of the corresponding lattices, de- 
noted FBZ(/ii) and FBZ(/i2)- It can then be shown that 
the pseudoinverse solution to the integral equation © is 
given by the inversion formula 

Sa(r) = / d 2 q / d 2 pK(r; q, p)0(q - p, p) , 

JFBZ(hi) JFBZ(ho) 



_ (5) 

where the kernel K is defined in Ref. |l2| . Several aspects 
of Eq. © are important to note. First, the transverse 
spatial resolution of reconstructed images is determined 
by the spatial frequency of sampling of the data function 
with respect to both source and detector coordinates. As 
a consequence, a large number of source-detector pairs is 
required to achieve the highest possible spatial resolution. 
It can be seen that when the source and detector lattices 
have equal spacing, the theoretical limit of transverse res- 
olution is given by the lattice spacing. When the source 
and detector lattice spacings are different, as is the case 
in the experiment reported here (where hi =0.5 cm and 
lri2 = 0.065 cm), the resolution of reconstructed images is 
controlled by the larger lattice spacing (lower spatial fre- 
quency). Second, the inverse problem in OT is evidently 
overdetermined. In addition, it is highly ill-posed. As a 
result, it can be said that large data sets allow for averag- 
ing the data function in such a way that the sensitivity to 
noise in the inverse problem is partially ameliorated. Fi- 
nally, numerical implementation of (0) requires replacing 
the integrals over d 2 q and d 2 p by sums over a finite set 
of wavevectors. In practice, we find that integration over 
d 2 q can be carried out with a step size Aq = 0.07 cm -1 
and 14,641 integration points while the integration over 
d 2 p requires a step size Ap — 1/2 Aq and 1,296 integra- 
tion points. Thus a total of 1.9 x 10 7 Fourier components 
of the data are used in the reconstruction. 

The first step in the reconstruction of tomographic im- 
ages is to measure the reference intensity 1$ for each 
source-detector pair. By fitting this data in the spatial 
frequency domain to (|3J) with a = ao we obtain the dif- 
fuse wavenumber ko = y^ao/D = 0.58 cm -1 and the 
extrapolation length £ = 0.7 cm. Note that these param- 
eters define the diffusion Green's function G in the slab 
geometry [l£| and that c^o and D cannot be separately de- 
termined from a continuous- wave measurement at a sin- 
gle wavelength. Next, the object to be imaged is placed in 
the sample chamber and the intensity / for each source- 
detector pair is measured. In Fig. 2 we show the recon- 
struction of a pair of black metal balls. The balls have a 
diameter of 8 mm and were suspended in the midplane of 
the sample chamber at a constant height with a separa- 
tion of 3.2 cm. Tomographic images were reconstructed 
with a 15 x 15 cm 2 field of view using 230 x 230 pixels per 
image with a separation between the slices of 0.26 cm. It 
can be seen in the central slice, which is equidistant from 
the source and detector planes, that the balls are well 
resolved. The shallower and deeper slices show that the 
balls remain well resolved but with a smaller diameter, as 
expected. Fig. 3 is a plot of 5a/ ao along the line passing 
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through the centers of both balls in the central slice. The 
distance between the peaks is 3.3 cm in close agreement 
with the measured separation of the balls. The FWHM 
of the peaks is 1.1 cm which slightly overestimates the 
diameter of the balls. The FWHM of the peaks in the 
depth direction is 1.5 cm (graph not shown). It is im- 
portant to note that the reconstructed contrast in 5a 
is not expected to be quantitative due to the possible 
breakdown of (0 in the interior of the strongly absorb- 
ing balls. Interestingly, however, the shape and volume 
of the spherical absorbers is recovered well. 
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Fig. 2. Reconstructions of 5a for the two-ball phantom 
plotted on a linear color scale. The distance of each slice from 



the plane of sources is indicated. All images are normalized 
to the maximum of the central slice. 

In conclusion, we have demonstrated the feasibility of 
analytic methods for image reconstruction in OT with 
large data sets. We are currently conducting further 
studies to assess the effects of absorption contrast on 
image resolution. In addition, the recent availability of 
CCDs with faster data acquisition will allow the collec- 
tion of data sets with greater numbers of sources, lead- 
ing to improvements in spatial resolution. We expect 
that with further technological advances, resolution con- 
sistent with the results of numerical simulations will 
be achieved. 
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Fig. 3. A one-dimensional profile of the reconstructed 
absorption along the line passing through the centers of the 
balls in the central slice. 
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